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ABSTRACT 

Context. Predicting the emerging X-ray spectra in several astrophysical objects is of great importance, in particular- when the 
observational data are compared with theoretical models. This requires developing numerical routines for the solution of the 
radiative transfer equation according to the expected physical conditions of the systems under study. 

Aims. We have developed an algorithm solving the radiative transfer equation in the Fokker-Planck approximation when both 
thermal and bulk Comptonization take place. The algorithm is essentially a relaxation method, where stable solutions are 
obtained when the system has reached its steady-state equilibrium. 

Methods. We obtained the solution of the radiative transfer equation in the two-dimensional domain defined by the photon 
energy E and optical depth of the system r using finite-differences for the partial derivatives, and imposing specific boundary 
conditions for the solutions. We treated the case of cylindrical accretion onto a magnetized neutron star. 

Results. We considered a blackbody seed spectrum of photons with exponential distribution across the accretion column and for 
an accretion where the velocity reaches its maximum at the stellar surface and at the top of the accretion column, respectively. 
In both cases higher values of the electron temperature and of the optical depth r produce flatter and harder spectra. Other 
parameters contributing to the spectral formation are the steepness of the vertical velocity profile, the albedo at the star 
surface, and the radius of the accretion column. The latter parameter modifies the emerging spectra in a specular way for the 
two assumed accretion profiles. 

Conclusions. The algorithm has been implemented in the xspec package for X-ray spectral fitting and is specifically dedicated 
to the physical framework of accretion at the polar cap of a neutron star with a high magnetic field (IS, 10 12 G). This latter case 
is expected to be typical of accreting systems such as X-ray pulsars and supergiant fast X-ray transients. 
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1. Introduction 

The solution of the radiative transfer equation (RTE) 
that describes the modification of a seed photon spectrum 
due to Comptonization in a plasma is a much debated 
mathematical problem. The equation in its full form is 
indeed integro-differential (Pomraning 1973) and allows 
for analytical solutions under some particular assump- 
tions, such as electron temperature T e — 0 (Titarchuk & 
Zannias 1998) or in the energy domain when the emerg- 
ing spectrum is a power law (Titarchuk & Lyubarskij 
1995). If the photon energy exchange for scattering is low 
(Au/u -c 1), it is possible to perform a Taylor expan- 
sion of the Comptonization operator around the photon 


initial energy, which transforms the RTE from integro- 
differential to purely differential (Rybicki & Lightman 
1979); this is known as the Fokker-Planck (FP) approx- 
imation. The necessary conditions to allow this mathe- 
matical approach are that the Compton-scattering process 
occurs below the Klein-Nishina regime, namely when the 
electron temperature kT c is subrelativistic (^ 100 keV), 
and that the optical depth of the Comptonization region 
is r 1. The regime of low temperature and high op- 
tical depth of the plasma indeed ensures that the spec- 
trum is almost isotropized, so that it is possible to use 
the Eddington approximation for the specific intensity 
J( v) = J(u) + 3V-P(i/), where J and F are the zero and 
first moment of the intensity field, respectively. 
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Moreover, if the plasma is not static but subject to 
dynamical (bulk) motion with velocity v(r), then another 
condition for using the FP approximation is that v(r) 
must be subrelativistic. When all the above restrictions 
are considered, one obtains by computing the first two 
moments of the RTE a main equation that describes the 
shape of the angle-averaged emerging intensity J(v) of 
the Comptonization spectrum. The general form of the 
RTE for the photon occupation number n(u) = J(u)/i / 3 
with subrelativistic electron temperature in the presence 
of bulk motion was first derived by Blandford & Payne 
(1981, hereafter BP81). Later, Titarchuk et al. (1997, here- 
after TMK97) showed that analytical solutions can be 
found if the velocity profile (assuming spherical symme- 
try) of the matter follows the free-fall law, vr oc R~ x / 2 . 
In this case, the equation can be written as L x n(x, r) + 
L T n(x,T ) = —s(x,t) and the solution can be obtained 
using the variable-separation method in the form 

OC 

n(x, r) = £ Cfc i? fc (r)iV fc (x), (1) 

k=l 

where c* and Rk{x) are the expansion coefficients and 
eigenfunctions of the space operator L T , respectively, 
while Nk(x) is the solution of the differential equation 

La;iVfc(x) — 'YfcfV).(x) = —s(x), (2) 

where 7* <x A| and A 2 is the k th -eigenvalue of the space 
operator. 

The Comptonization spectrum is mostly dominated 
bv the first eigenvalue (see Fig. 3 in TMK97), while the 
terms Njt (x) with k > 2 represent the fraction of pho- 
tons that leave the medium without appreciable energy- 
exchange. Starting from the results of TMK97, Farinelli 
et al. (2C08, hereafter F08) developed a model (comptb) 
for the X-ray spectral fitting package xspec, which com- 
putes the emerging spectrum b3' means of a numerical 
convolution of the Green’s function of the energy operator 
with a blackbody (BB)-like seed spectrum. The model has 
been successfully applied to a sample of low-mass X-ray 
binaries hosting a neutron star (NS). The method of the 
variable separation has also been adopted by Becker & 
Wolff (2C07, hereafter BW07), to find analytical solutions 
of the RTE in the case of cylindrical accretion onto the 
polar cap of a magnetized NS. The starting equation in 
BW07 is formally the same as in BP81: the most signif- 
icant difference is that the Thomson cross-section is re- 
placed by an angle-averaged cross-section that takes into 
account the presence of the magnetic field ( B ~ 10 12 
G). Following the results of Lyubarskii & Sunyaev (1982), 
BW07 assumed a velocity profile v(t) oc — r, which al- 
lowed the RTE to be separable in energy and space. Note 
that the assumed velocity profile implies that the matter 
flow stagnates at the stellar surface, which is at odds with 
the solution of TMK97, where the matter velocity is in- 
creasing ;owards the central object, which can be either 
a NS or a black hole (BH). When the velocity profile is 


not free fall-like (TMK97, F08), or oc r (BW07), the vari- 
able separation method can no longer be applied and the 
solution of the RTE can be obtained only with numerical 
methods. 

We report a numerical algorithm that allows the solu- 
tion of the RTE in the FP regime using finite-differences 
for any desired velocity profile and seed photon spatial and 
energy distribution. We apply in particular it to cylin- 
drical accretion towards the polar caps of a magnetized 
NS, following the approach of BW07. The algorithm es- 
sentially uses a relaxation method, therefore it finds the 
asymptotic (stationary) solution of the RTE for a given 
initial value (at time t = 0) condition. Our work is struc- 
tured as follows: in Section 2 we describe the kernel of the 
algorithm for generic two- variable elliptic partial differen- 
tial equations; in Section 3 we formulate the problem for 
the general RTE and appropriate boundary conditions; in 
Section 4 we consider the more specific case of a system 
configuration with azimuthal symmetry, typical of cylin- 
drical accretion; in Section 5 we show the emerging spec- 
tra obtained for different sets of the theoretical parameter 
space; finally, in Section 6 we briefly discuss possible as- 
trophysical consequences and implementations (e.g., for 
xspec) derived from the application of the algorithm. 


2. The general elliptic partial differential 
equations 

The algorithm we report is essentially based on the relax- 
ation methods, which allow one to find the solution of a 
boundary elliptical problem. The differential equation has 
to be written by finite differences. Once the sparse matrix 
is defined, it can be split into layers over which an itera- 
tion process is applied until a solution is found (Press et al. 
1992). The general form of a linear second-order elliptic 
equation with vanishing mixed derivatives and a source 
term can be written as 

^ + Q(x,y)^ + R{x,y)u + w{x,y)^ 

+Z(x,y)^ = ^-S{x,y). ( 3 ) 

We define a three-dimensional grid of discrete points for 
the variables x, y and f ; 

Xi = xo 4- ih x , i = 0, 

Vj = I/O + jky , j = 0, 1, . . . , Nyt 

t m — to+ mh t , m = 1, 2, . . . , M, (4) 

where h x ,h v ,ht are the grid spacing. The function 
u(x, y, t) is evaluated at any point of the grid, so we write 
it as ui m . We write the first and second derivatives over 
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the variables using finite differences: 


the system (9), which gives 


du 


d 2 u 


2ur + u j r 

dx 

h x 

dx 2 


hi 

f 

du 

4 +i ’ m - u r 

d 2 u 

_4 +1 

,m 

2ui' m + 4 

dy 

hy 

W 



hi 

du 

u j,m _ u J,m- 1 





dt 

h t 






A x u m + A y u m = 


u — u 


m— 1 


h t 


- S 


+ht A x A y ( u m - u 17 *- 1 ) . 


( 10 ) 


The third term on the right-hand side of equation (10) 
represents the residual error in the numerical solution. As 
(5) a first step, we collect the terms with the same index in 
both equations (9) and obtain 
Substituting the above definitions into equation (3) and , ^ ^ 

collecting terms, we obtain f A x - — J w m_1 / 2 — - ( A y + — J u m ~ 1 — S, 

+ ^ + ^- 1 ' m + + ^ +1 ' m / ‘ x \ , 1—1/2 

' ‘ ’ ,+1 * (A v -pU m = A . (11) 


ht 


-si, 


where 

n j _ 'P[xi,y j ) 

1 hi ’ 


fh 

Both equations are defined inside a 2D (x, j/)-domain, 
(®)with boundary conditions defined according to the spe- 
cific problem under consideration. First, for any m and 
j values, we must impose the boundary condition on the 
left-hand side of the x-domain ( i = 0) for the function 


u\ 


j,m- 1/2 


2V(xi,y j ) Q(xi,y J ) 


<£ = 


fi 


hi 

hx 

V{xi,y j ) , 

Q{xi,y j ) 

*1 

h x ’ 

W(xi,V) 


h 2 ’ 
n y 


2W(x u yi) 

Z{x,,y j ) 

K 

hy 

W(xi,yi) 

Z{xuyj) 






= So. 


(12) 


while the source term Sq is defined at the beginning. Thus, 
for i = 0, equation (6) can be written as 


^r'/ 2 =L&r- l,2 +Ki 

where 


(13) 


Li 


4 


y _ _L 
°o ht 


Qj,m - 1 

fri - 

K ° ~ y _ A 

°o h t 


K 


hy 


S{ = S(x t ,y j ). 


(7) 


sr - 1 = - (4 +4+ ^ + i) u ^- 1 - 


(14) 


, ., , ... c , with the coefficients determined m equation (7). For i = 1, 

The operators over the x and y vanaoles are then defined . , v ' ’ 

usmg equation (13) we obtain 


as 


A a, — a{ + bj 

Ay = 4 + 4 +fi. 


(8) 


(a? Ll + 4 ^ju{ m ~ 1/2 +44 


ml/2 


The solution procedure consists of dividing equation (6) which can be written M 
into two equations. The first gives us the solution for an 
intermediate m — 1/2 layer, while the second provides the 
solution for the m-layer. As starting point we need to es- 
tablish an initial guess for the function at m — 1. The w ^ icre 
system of equations to be solved is thus 


,.m- 1/2 _ 1 

A x u m_1>/2 + \ y u m ~ x = — S, 


L[ = - 


= s{ m ~ x -4ki, 

IS 

r- 1 / 2 = L{ u L m - i/2 + k\, 

4 


(15) 

(16) 


i . j frj 

1 a l A 0 /i 7 \ 

4~r t ' 1 a {Lo + 4-t 


A y(u m — U m ‘) = 


m u m- 1/2 

ht ’ 


(9) 


a{L J 0 + 1 

Iterating the process, we obtain the general form 

nj-.*— i/3 =ti4£r i / 2 +Ki, 


(18) 


C? 


in which we have temporarily dropped the indices i,j. where 
From the system we notice that the intermediate layer 

m — 1/2 is needed only to build up the solution at the Lj = — pr— - . 1 

subsequent layer in m. The numerical accuracy of the so- °i^»- 1 % — hi 

lution can be estimated by combining both equations in 


ki - 


sr~ x - ajkj_ 1 

4 lLi + 4 -k' 


( 19 ) 
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At the right boundary of the x-domain (i — N x ), we im- 
pose the second boundary condition 

(»> 

Nov/, using equation (18) we can thus build up the solution 
over the x-variable iteratively as 

j,m-l/2 fj J fcj 
U N X - 1 ''JVj-IWj + n N x -l< 

j,m- 1/2 fj J,m—l/2 f>j 
u N x —2 ~ lj N x ~2 u N x -l ~ n N x -2’ 

uj m - 1/2 = Lj«r^ + 4 (21) 

Therefore, the construction of the solution is obtained in 
two steps: a bottom-up process which allows one to build 
the coefficients L\ and K\ (Eq. [19]) starting from the left 
boundary condition on u 3 0 (Eq. [12]), followed by a top- 
down procedure determined by the right boundary condi- 
tion v? Nx (Eq. [20]). 

Once the solution over the x- variable for the m — 1/2 
layer is obtained for any j (the index of the y variable), we 
then seek the solution of the second equation in the system 
(9) by following the same procedure described above with 
initial boundary condition at j = 0 

u°' m = L?u)' m + Kf, (22) 

and, similarly to equation (19) 


dJZr + e? + £ 


Sl' m -d\K{- 


v'here 


/ i \ j ,m- 

sr = (4 + 4 + /1 + 1) 4 m - 1 - \ 


j>m— 1/2 


depends on the solutions u J i ' nt ~ 1 ^ 2 and obtained in 

the layers m— 1/2 and m — 1. 

As required for the procedure over the X- variable, the coef- 
ficients L\ and K\, built from J? t and K°, are determined 
by the left boundary condition (j = 0) for the function 
u/, and the solution for any j is determined by the right 
boundary condition uf" = g^’ y : 

= L^~ l u" v ‘ m + K? v ~\ 
u Ny-2,m = L N v -2 u N v -l,m + ^N v - 2, 

u°’ m = + K°. (25) 

After constructing the solution over the x and y variable, 
the solution of the top layer m becomes the initial function 
of the bottom layer related to iteration m + 1 according 
to the scheme 

m = 1 — > (u°, u 1 / 2 , u 1 ), 
m= 2 — > {u l ,u^ 2 ,u 2 ), 


and u 3 / 


obtained in 


It is also worth mentioning that at the first iteration m 1 
an initial guess function u{'° must be assumed, which is 
needed to find the solutions u\ }1 ^ and uj’ 1 in the system 
(9). The loop over m stops when the same convergence 
criterion is satisfied, which physically means that the so- 
lution has "relaxed"’ to its stationary value. One possible 
criterion could be 1 — £ < \u\ m ~ l \/\u\ m ~ 1 ^ < 1 + £ and 
1 — £ < [u| m-1/,2 |/|u| m < 1 + £, where £ is a user-defined 
numerical tolerance. In the next section we will show in- 
stead another convergence criterion we have chosen for 
stopping the iteration procedure, for the particular case 
of the RTE. 

3. Application to the radiative transfer equation 
and boundary conditions 

The general form of the RTE in the presence of subrela- 
tivistic bulk motion for a plasma with constant tempera- 
ture T e is given by (see Eq.[18] in BP81) 


■+V- Vn = V- 


,3n e a(u) 


\ 1 /_ ,,, dn 

.) + 3<V-V)^ 


1 d [n e tr(r') 4 / dn\ .. 


where n(v, r) is the zero-moment occupation number of 
the radiation field intensity, V is the plasma bulk veloc- 
ity vector, a(v) is the electron scattering cross-section, 
ne(r) is the electron density and j{ u, r) is the source term. 
Because the spectral formation is determined by the op- 
tical depth r of the system, we use the latter quantity as 
the actual space variable. The solution of equation (27) is 
fulfilled by imposing the boundary condition at the sur- 
face defined bv r = 0, (which represents the starting point 
of the integration domain) for the spectral flux, which is 
given by 

¥{v,r) = -v 3 \( 1 , , Vn) + . (28) 

Under particular symmetries of the system configuration 
(e.g., cylindrical or spherical), the problem becomes one- 
dimensional. For constant electron temperature T e it is 
also more convenient to use the adimensional variable 
x = hi>/kT e \ moreover, when performing numerical inte- 
gration using finite-difference methods, we use a logarith- 
mic binning of the energy through the additional change 
of variable x — ► e g . Under these assumptions, equation 
(28) becomes 

F ^‘-[lTr + l V {^ - J )]' < 29 > 

where J = n x 3 is the specific intensity. 

At the inner boundary we impose the condition 


m = M -> (u M “ 1 ,u M-1 / 2 ,u M ). 
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where A is the albedo at the surface. A fully absorptive 
surface (A = 0) is appropriate for a BII, while 0 < A < 1 
accounts, e.g., for a NS atmosphere. However, the in- 
ner boundary condition (30) depends on energy as well 
as space (see Eqs. [28] and [29]). For mixed boundary 
value problems, no analytical solutions are possible (see 
Appendix E in TMK97) and numerical methods prove 
to be unstable. However, in the energy range where the 
spectrum is a powerlaw J(x,t) = R(r)x~ a , equation (30) 
becomes 

-f + A(a + 3 = (3D 

where ft, is the bulk velocity at the inner radius (r = 0), 
and here the problem is reduced to a standard boundary 
condition over the space variable r. Writing the derivative 
in terms of finite-difference, equation (31) then becomes 


is above the major contribution of the BB component and 
below the expected high-energv cut-off value. 

Once a m is estimated, it is inserted into equation (34), 
which accordingly represents the boundary condition at 
t = 0 for the iteration m- 1-1 (see Eq. [22]). We then com- 
puted the new index Q m +i for and again inserted 

it into equation (34) at iteration m + 2, and so on. The 
routine is stopped when a m and a m +i differ less than 
10 -5 provided that the condition holds for a sufficiently 
high number of iterations (> 100). Note that the same 
criterion is adopted also if ft = 0, even if then of course 
L° remains constant across the iteration. We have also 
verified that this criterion automatically also satisfies the 
convergence of the norms |u| TO_1 , |u| m-1 / 2 , and |u| m . 

4. Cylindrical accretion onto a magnetized 
neutron star 


u] — u? . . 0 3/1 — 0 

— ir +A<a+3) “‘““ 2 (it+K 

which can be written after collecting terms as 

1 


«? = 


u 


l + h T [0 o (a + 3)+G{A)} * 


(32) 


(33) 


where G(A) = 3/2(1 — A)/{\ + A). We then set our prob- 
lem as follows: first, because the solution u\ of the RTE 
physically represents a specific intensity, it must by defi- 
nition be equal to zero in the limits E -> 0 and E -4 oo, 
therefore we set v? 0 = u 3 N =0 (see Eqs. [12] and [20]). For 
the behavior of the function uj for r = 0 (J = 0), equation 
(33) immediately allows us to define (see Eq. [22]) 



1 

1 4- ft [ft (a 4 3) + G(A)] 


K? = 0. (34) 


N 

As outer boundary condition over r, we impose that u i v = 
0, which means that the specific intensity goes to zero for 

X ^ 'T max . 

We emphasize that the condition u\ > 0 for any 
value implies a specific restriction in the choice of the step 
size h T , which ensures that L® > 0 (as ft < 0). More 
specifically, we imposed the condition on h T such that the 
number of steps over r be N~ = r max /h T > 10. 


3.1. The iteration procedure 

As already mentioned in Section 2, it is necessary to choose 
a convergence criterion for stopping the iteration over the 
m variable. We proceeded in the following way: at each 
iteration m, we computed the spectral index a m of the 
solution (corresponding to r = 0) in a given energy 
range E m \ n — E max . To minimize bias or wrong estimate 
of a, 7 1 , the definition of the energy interval for the com- 
putation of the spectral slope must be chosen carefully. If 
the seed photon spectrum is a BB with temperature kTbb, 
a reasonable choice can be the assumption Emin ~ 7kTbb 
and E max « 20A;Tbb, respectively, given that this interval 


We applied our algorithm to solve the RTE for accretion 
towards the polar cap of a magnetized NS, whose mathe- 
matical formalism was developed by BW07 in the frame- 
work of the spectral formation of accretion-powered X-ray 
pulsars. The relatively strong magnetic field (B 10 12 G) 
of the NS is expected to channel the accretion flow to- 
wards the poiar caps, and for low values of the altitude 
above the star surface, the problem can be treated in a 
axis-symmetric approximation where the space variable is 
defined by the vertical coordinate Z. The magnetic field 
moreover forces the medium to become birefringent as 
the effect of vacuum polarization, and birefringence en- 
tails the formation of two linearly polarized modes (or- 
dinary and extraordinary) of the photons, each having a 
characteristic scattering cross-section. For ordinary mode 
photons with energy below the first cyclotron harmonic 
at E c « 11.57 Si 2 keV (where B i2 = B( 10 12 G), BW07 
defined angle-averaged cross-sections parallel and perpen- 
dicular to the lines of the magnetic field as crj j = 10 -3 <tt 
and o± = ax, respectively, where or is the Thomson scat- 
tering cross-section. This is indeed the only approximation 
that allows to treat the problem analytically or numeri- 
cally. We note that Ferrigno et al. (2009), starting from the 
analytical solutions reported in BW07, developed a model 
that was later almost successfully tested on the accret- 
ing pulsar 4U 0115+63. Their model is based essentially 
on the convolution of the column-integrated Green’s func- 
tion of the thermal plus bulk scattering operator with a 
given seed photon distribution. The basic assumption of 
this derivation is that the velocity profile of the accreting 
matter is assumed to be v(t) ex — r, which allows one to 
find analytical solutions through the variable separation 
method (Eqs. [36] and [37] in BW07). The numerical al- 
gorithm we developed directlv solves the RTE, without 
the need of this prescription for the dynamical configura- 
tion of the accreting matter field, and we included some 
modifications with respect to the approach of BW07 and 
Ferrigno et al. (2009). 

First, following TMK97, we include in equation (27) a 
second term in the thermal Comptonization operator that 
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accounts for the contribution of the bulk motion veloc- 
ity of electrons in addition to their thermal (Maxwellian) 
component. With this prescription in mind, equation (27) 
becomes 


To solve equation (37), it is necessary to define the behav- 
ior of the velocity profile &{r). We considered two possi- 
bilities: in the first one, we assumed a general form 


P{Z) = 

where the normalization constant 


-*) 


(42) 

defined = 


l Q n v dn dv e dn d ( 1 dn\ where the normalization constant is denned .& = 

cm ~ S{ - e,Z) = ~~cdZ + dzTcm + dZ (^Tdz) Po{Z 0 /Z ,)", and A is the terminal velocity at the alti- 

V 11 ' tude Z 0 . 

(35) The continuity equation for the system here considered 
gives the electron number density 


n 


+ 


n e a 1 d 
m e c 2 e 2 dc 


e 4 ( n + (kT c + m e v 2 /S) 


dn 

m 


)]' 


where c = hu, a = 10 1 <tt, while i esc is the photon mean 
escape timescale (see Eq. [17] in BW07) 


teBC — 


n e a±r% 


1 dJ S(q, t) 


n e cr\\cH m H 


1 + 


m e v(r) 


21 


+ 


+ 


e q -35- 


He 2 


.21 


3fcT e J dq 2 

dJ_ 
dq 


1 d 2 J v(t) dJ 

+ 3 Urn 2 Hem' 


3 kT e (e g — 3 + 5) — m e v(r) 2 


3JtT e 


where we have defined the quantities 

a kT e 


and 


H - 


<5 = 


a || m e c 2 ’ 
1 dP(T) 


c = 


m 


3kT e (e q — 3 + 5) — m e u(r) 2 
3 kTe 
£ 2 v{t) 2 


Q(r) = 

TI(t) = e 9 - 3 6 


He 2 ’ 




Z(t) = - 
S{q,r) = 


v(j) 
He ’ 
S{q, t) 
H 


TLq — 


M 


7rm p | J 0(Z’)|ci?o ’ 


(43) 


(36) 


Now, using the relation dr — n e oudZ and the logarithmic 
binning of the adimensional energy x = hv/kT e , equation 
(35) becomes 


(37) 


(38) 


(39) 


.where m = M/Me is the mass accretion rate in Eddington 
units and Ro is the radius of the accretion column. 

We then define the adimensional quantities z and ro 
through the change of variables Z -+ 7?sg mz and Ro —> 
i?s©mro, where m = M/Mg, while Ms© and f?s© are 
the Sun mass and Schvarzschild radius, respectively. The 
effective vertical optical depth of the accretion column is 
then given by 

r m ( zV+1 - 

r( z ) = j^dZ' = C— 2 — — , (44) 

where C = 2.2 x 10^ 3 , and zo is the vertical coordinate at 
the NS surface. 

Inverting relation (44), we also define the velocity profile 
of the accreting matter as a function of the optical depth 
t instead of the space variable z 


0(t) 




0+1 + -^ r o(l + i?)t ) ” rl 


Cm 


(45) 


3 H dr ’ 

where 0 It) — v(t)/c, while the dimensionless parameter 
/ is given by (see Eq. [26] in BW07) 

15.8 r 0 


As a second possibility, following BW07, we considered 
the velocity profile 


P( T ) = 


(46) 


(40) 


Equation (37) is given in the general form (3) and for this 
particular case, we have 




( 41 ) 


where $ = 0.67 £/z 0 (see Eq. [32] in BW07). 

Given that in our model the optical depth r represents 
one of the free parameters, once it is provided in input 
together with adimensional accretion column radius ro, 
the accretion rate m must be first computed either from 
equation (44), if 0{r) is defined as in equation (45), or 
from equation (28) in BW07 if 0(r) belongs to equation 
(46). This step is necessary to determine the £ parameter 
(Eq. [40]), and requires fixing the maximum altitude of the 
accretion column z max . We assumed z max ~ 2zq, and all 
emerging spectra (see next section) were computed with 
this choice. 

5. Results 

In this section we report some examples of the theoreti- 
cal spectra obtained by the numerical solution of equation 
(37) for different sets of the physical quantities that define 
the system. We consider a BB seed photon spectrum at 
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given temperature fcTbb with exponential spatial distribu- 
tion across the vertical direction, according to 


S(x, r) = C n e~ 


kT?x 3 

gfcT'e/feT'bb * — . 


(47) 


with the normalization constant defined as C n = 
where 7?km and £>io are the BB emitting area in 
kilometers and the source distance in units of 10 kpc : re- 
spectively. The spectra were computed using the velocity 
profiles defined in equations (45) and (46), respectively. 
The common parameters for both cases are consequently 
the BB temperature fcTbb, the electron temperature kT e , 
the optical depth r, the albedo at the inner surface A 
and the radius of the accreting column ro. On the other 
hand, for /3(r) belonging to equation (45), additional pa- 
rameters are the index 77 and the terminal velocity at the 
star surface 0o- We first present the results for this second 
physical case. 

In Fig. 1 we show the emerging spectra for different 
values of the electron temperature kT e and two terminal 
velocities 0o = 0.1 and 0a = 0.64. As expected, both times 
higher values of kT e produce flatter spectra and push the 
cut-off energy E c to higher energies; on the other hand, the 
bulk contribution as a second channel of Comptonization 
depends on the value of kT e . The two extreme temperature 
values reported here, kT e — 5 keV and kT e = 50 keV, 
are particularly instructive: for low electron temperatures 
the spectrum changes from BB-like when 0q = 0.1 to a 
cut-off power law with E c 30 keV when 0a — 0.64, 
while the spectral change is much less enhanced for a hot 
plasma. These can be considered as typical examples of 
bulk-dominated and thermal-dominated Comptonization 
spectra, respectively. 

Together with the electron temperature, the optical 
depth r is an important parameter that plays a key role 
in determining the spectral slope and cut-off energy, as 
clearly shown in Fig. 2. We note that in Fig. 1 and Fig. 2 
the index of the velocity profile was chosen to be rj = 0.5, 
typical of accretion onto a compact object where gravity 
and radiation pressure are the only force terms that de- 
termine the dynamical configuration. Here, the terminal 
value of the matter velocity 0q depends on the ratio of 
the radiative and gravitational forces, provided the con- 
dition |F r |/|F g | £ 1 is satisfied. This relatively simple 
approach is valid for low values of optical depth r, while 
when r > 1 radiative transfer becomes important and the 
problem requires in principle a more accurate radiation- 
hydrodynamics treatment. 

It is outside the scope of this paper to compute the 
exact velocity profile for accreting matter under the pres- 
ence of a strong radiation field in a high optical depth 
environment. We merely introduced a simple parametriza- 
tion for modifying the velocity field by changing the index 
77 , with Lie results shown in Fig. 3, for two different val- 
ues of the electron temperature kT e . As Fig. 3 shows, the 
lower the value of 77 , the harder the spectrum: this behav- 
ior can be explained in a quantitative and a qualitative 


way. Indeed, as 77 increases, the velocity profile 0(z) be- 
comes sharper, and for a fixed terminal velocity 0o, elec- 
tron temperature kT e , and optical depth r, while photons 
diffuse through the bounded medium, on average the en- 
ergy of the electrons (caused by their Maxwellian plus 
bulk motion) will be lower, and consequently the net en- 
ergy gain of the photons due to inverse Compton will be 
less. From the mathematical point of view, it is worth 
mentioning that Mastichiadis & Kylafis (1992, hereafter 
MK92) reported the analytical solution of the RTE in 
the Fokker-Planck approximation with the variable sep- 
aration method for spherical accretion without magnetic 
field in the limit T e = 0 . Assuming a general velocity pro- 
file 0 X oc 7-~ t) , the authors showed that the spectral index 
of the fc th -Comptonization order emerging spectrum yields 
Ok = 3 + 3Ak/(2 — 77 ) (see Eq. [1]), where Ak is the fc th - 
eigenvalue of the space operator. Using equation (B12) of 
MK92, it follows immediately that as 77 increases, the spec- 
tral index Ok increases as well. This mathematical result 
in terms of spectral formation can be considered as general 
in the framework of the FP treatment, and is accordingly 
qualitatively meaningful for our research. 

We also emphasize that analytical solutions for 77 0.5 

have been possible for MK92 only because of the condition 
T e = 0, which drops the thermal Comptonization operator 
in the RTE, while when T e > 0 this is possible only for 
77 = 0.5 (TMK97, F08). 

In Fig. 4 we show results for different terminal bulk ve- 
locities 0o for two electron temperature values. The figure 
can be considered an extension and completion of Fig. 1 
because more values of 0o are shown to better appreciate 
the induced changes in the emerging spectra. 

The spectral modifications as a result of different val- 
ues of the albedo A at the inner surface are instead shown 
in Fig. 5, where we explored full absorption (A = 0) and 
full reflection (A = 1), together with other intermediate 
values. In the framework of a physical link to astrophysical 
objects it would be natural to associate a BH to the condi- 
tion A = 0 and a NS to the condition A = 1 , respectively 
(Titarchuk & Fiorito 2004; Farinelli & Titarchuk 2011), 
even though this latter assumption may be considered 
an oversimplification of the problem. A most realistic ap- 
proach would consist indeed in an energy-dependent treat- 
ment of the albedo, a problem that could be faced only 
with Montecarlo simulations, with the additional compli- 
cations arising from a detailed treatment of the star pho- 
tosphere (surface) properties. 

For our unavoidably simplified assumptions, the net 
effect of increasing values of A is a progressive flattening of 
the emerging spectra. This is physically explained because 
when A > 0, a fraction of photons (which becomes 100% 
when A=l) suffers on average more scattering with respect 
to A = 0. Qualitatively, the spectral modification leads in 
the same direction as an enhanced optical depth of the 
system. 

The last parameter that strongly influences the spec- 
tral formation is the radius of the accretion column ro, 
whose effects are shown in Fig. 6 . Indeed, following the 
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BW07 prescription, the mean escape time for photons 
using the diffusion approximation tesc oc (see Eq. 
[36]). On the other hand, both the bulk and thermal 
Comptonization parameters (ybuik and y t h, respectively) 
are related to the mean number of scatterings that pho- 
tons experience in the medium via 

{/bulk » iV^ Ik Cbuik, (48) 

J/th ^ f^avCth) 

where N^ lk , Cbuiki Nlv an 4 Cth are the averaged number 
of scatterings and the fraction energy gain per scattering 
for bulk and thermal Comptonization, respectively. Both 
N^v lk and IVj* are of course also proportional to i esc (see 
Eqs.[94]-[97] in BW07). 

Evidently therefore, for fixed velocity profile parame- 
ters srf and r] (see Eq. [42]), once the optical depth r is de- 
fined (see Eq. [44]), to keep its value constant for increas- 
ing ro (as reported in Fig. 6), the accretion rate rh must 
also increase in a way to keep the ratio m/r q constant. 
Combining equation (36) and (43) yields t e 5C oc rh, which 
in turn leads to an enhancement of the Comptonization 
parameters j/buik and j/ t h in equation (49) with a harden- 
ing of the spectral shape. 

Considering now the velocity profile defined in equa- 
tion (46), we see that the results are qualitatively the same 
as in equation (45) as far as the spectral modifications in- 
duced by variations of kT e are concerned (Fig. 7), r (Fig. 
8) and A (Fig. 9), respectively. But there are opposite ef- 
fects that are induced in the emerging spectra by different 
values of the accretion column radius ro for the velocity 
profile here considered. Indeed, using equation (40) and 
the definition of r in equation (28) of BW07, which al- 
lows us to express the accreting matter velocity in terms 
of the z-coordinate, in spite of the optical depth r, it 
is straightforward to see that fi(z ) oc r 0 '^ 2 . In particu- 
lar, if zq = 2.42 and z max = 2zo we have /3 ma *= 0.60 for 
r 0 = 0.1, jd max =0.38 for r 0 = 0.25, /3 max =0.27 for r 0 = 0.5 
and /9 ma x=0.2 for ro = 1, respectively. Note that because 
equation (46) describes matter that stagnates at the star 
surface, here /9 max represents the velocity at the accretion 
column altitude z max . In other words, while using equa- 
tion (45), the choice of ro does not modify the velocity 
field of the accreting matter, which is only determined by 
the choice of ft and r?, for (46) as ro increases the bulk 
contribution to the spectral formation becomes less impor- 
tant, and this drop is not compensated for by the increase 
of the photon mean escape time f esc , which, as explained 
above, would instead contribute to spectral hardening. 

6. Implementation in the xspec package 

Our model will be publically available and distributed as 
a contributed model to the official XSPEC 1 web page. 

In Table 1 we report a summary of the free parameters of 
the model with their physical meaning. The code is written 

1 http://heasarc.nasa.gov/xanadu/xspec/newmodels.html 


Table 1. Parameter description of the xspec model COMpmag. 


Parameter 

Units 

Description 

kTbb 

(keV) 

Seed photon blackbody temperature 

kT e 

(keV) 

Electron temperature 

T 


Optical depth of the accretion column 

n 


Index of the velocity profile (Eq. [45]) 

ft 


Terminal velocity at the NS surface 
(Eq. [45]) 

ro 


Radius of the accretion column in 
units of the NS Schwarzschild radius 

A 


Albedo at the NS surface 

Flag 


= 1, P{t) from equation (45) 
= 2, 8{t) from equation (46) 

Norm 


RlJDt o 


using C-language, and can be easily installed following the 
standard procedure reported in the official xspec manual 
and in the brief cookbook, which will be delivered together 
with the source code. As a general concern for users, we 
point out that usually the emerging spectrum obtained 
from the Comptonization of a seed photon population with 
any given energy distribution S(E) can be presented as 
the sum of the seed spectrum and its convolution with 
the scattering Green’s function G(E,Eo) of the electron 
plasma, each with their relative weight, according to the 
general formalism 

F(E) = -^~[S{E) + Ax S(E) * G(E , E 0 )), (49) 

where C n is a normalization constant. The ratio A/(A+ 1) 
is the Comptonization fraction, and its value qualitatively 
determines the contribution to the total spectrum of the 
Comptonized photons. The value of A, here not to be con- 
fused with the albedo, may depend on several geomet- 
rical and physical factors, such as the spatial seed pho- 
ton distribution inside the system configuration (see Fig. 
4 in TMK97). The lower the value of A, the more en- 
hanced the direct seed photon spectrum S(E). Examples 
of xspec models that use the definition in equation (49) 
are BMC (TMK97) and comptb (F08). Either model, how- 
ever, does not solve the full RTE including the photon 
spatial diffusion and distribution, the latter of which is an 
unknown quantity that is phenomenologically described 
through the continuum parameter log(A). On the other 
hand, it is not possible to change the value of log(A) ar- 
bitrarily in our present model, i.e., according to the ob- 
served spectra. Its value is implicitly determined once the 
seed photon spatial distribution is fixed. 

We presented the results of simulated spectra assum- 
ing an exponential distribution over r for S(E), which 
was assumed to be a BB; then, the transition from the 
low-energy part of the spectrum (the Rayleigh regime for 
E 3fcTbb) to the high-energy (Comptonized) powerlaw 
shape is almost smooth, which corresponds approximately 
to A » 1 in equation (49). Other seed photons spatial 
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distributions can produce a different onset between the 
BB peak and the powerlaw-like regime. In general, for 
observed spectra where a direct and enhanced BB-like 
component is required by the fit, our claim is to model 
the source continuum with modelization of the type BB — 
compmag by preferably keeping equal to each other the 
temperatures of the direct and Comptonized BB compo- 
nent. 

7. Conclusions 

We have developed a numerical code for solving elliptic 
partial differential equations based on a relaxation method 
with finite differences. In particular, we reported a spe- 
cific application of the algorithm to the radiative transfer 
equation in the Fokker-Planck approximation, which is of 
particular interest for high-energy astrophysical applica- 
tions. We considered cylindrical accretion onto the polar 
cap of a magnetized neutron star, using the mathemati- 
cal formulation of the problem reported in Becker & Wolff 
(2007) with some modifications. Specifically, we included 
the second-order bulk Comptonization term in the scatter- 
ing operator and we considered different velocity profiles 
for the accreting matter. 

The code for the computation of the emerging spec- 
tra in this configuration has been written with the aim 
to implement it in the xspec package and will be deliv- 
ered to the scientific community. Because angle-averaged 
cross-sections caused by the magnetic field were included, 
the model is suitable to be applied to the observed spectra 
of sources where most of spectral formation is claimed to 
form close to a NS with high magnetic field ( B 10 12 G), 
such as accreting X-ray pulsars and supergiant fast X-rav 
transients. Of course there are some unavoidable simplifi- 
cations in the model, such as the assumption of constant 
electron temperature of the accretion column. We note, 
however, that the isothermal condition is typical of any 
Comptonization model implemented in XSPEC because it 
is fundamental for users to reach a compromise between 
the accuracy of the physical treatment and the computa- 
tional speed. More accurate theoretical investigations of 
the accretion processes can be performed only with MHD 
simulations performed through PC cluster-computing re- 
sources, which are beyond the scope of the present work. 

As in many theoretical models, the number of available 
free parameters is higher than the number of the observ- 
able ones. Therefore a correct working approach is to keep 
some parameters fixed during the X-ray spectral fitting 
procedure to avoid degeneracy. 

If this model is used in any publication, we kindly ask 
the authors to cite this paper in the reference list. 

Acknowledgements. We acknowledge financial contribution 
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Fig. 1. Emerging spectra obtained from the solution of equation (37) for different values of the electron temperature kT e and 
velocity profile defined in equation (45). In both cases the fixed parameters are fcTbb = 1 keV, r = 0.2, rj = 0.5, no = 0.25, 
A = 1. Lift panel: /3o = 0.1. Right panel: fin = 0.64. 




Fig. 2. Same as Fig. 1 but for different values of the optical depth r. Fixed parameters are fcTbb = 1 keV, jj = 0.5, fin = 0.64, 
ro = 0.25, A = 1. Left panel: kT e = 5 keV. Right panel: kT e — 15 keV- 




Fig. 3. Same as Fig. 1 but for different values of the index of the velocity profile. Fixed parameters are fcTbb = 1 keV, r = 0.2 
fio = 0.64, ro = 0.25, A 1 . Left panel: kT e = 5 keV. Right panel: fcT 0 — 15 keV. 
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En&rgy (ksV) 



Fig. 4. Same as Fig. 1 but for different values of the inner velocity /Jo- Fixed parameters are fcTbb = 1 keV, r = 0.2, ij = 0.5, 
r 0 — 0.25, A = 1. Left panel: kT e = 5 keV. Right panel: kT e = 15 keV. 




Fig. 5. Same as Fig. 1 but for different values of the albedo A, with the velocity profile of equation (45). Fixed parameters are 
fcTbb - 1 keV, t = 0.4, rj -= 0.5, do = 0.64, ro — 0.25. Left panel: kT e — 5 keV, Right panel: kT e = 15 keV. 




Fig. 6 . Same as Fig. 1 but for different values of the accretion column radius ro, with the velocity profile of equation (45). Fixed 
parameters are fcTbb - 1 keV, r = 0.2, 77 = 0.5, do = 0.64, A = 1. Left panel: kT e = 5 keV. Right panel: kT e = 15 keV. 
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Fig. 7. Emerging spectra obtained from the solution of equation (37) for different values of the electron temperature kT e , with 
the velocity profile of equation (46). In both cases the fixed parameters are feTbh = 1 keV, ro = 0.25, A = 1. Left panel: r = 0.2. 
Right panel: r = 0.4. 




Fig. 8. Same as Fig. 1 but for different values of the optical depth r, with the velocity profile of equation (46). Fixed parameters 
are fcTbb = 1 keV, ro = 0.25, and A = 1. Left panel: kT e = 5 keV. Right panel: kT e — 15 keV. 




Fig. 9. Same as Fig. 1 but for different values of the albedo A, with the velocity profile of equation (46). Fixed parameters are 
kTbb - 1 keV, r = 0.4, and ro - 0.25. Left panel: kT e = 5 keV. Right panel: kT e = 15 keV. 
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Fig. 10. Same as Fig. 1 but for different values of the accretion column radius rg, with the velocity profile of equation (46). 
Fixed parameters are k7'bb ~ 1 keV, r = 0.2, and A = 1. Left panel: kT e = 5 keV. Right panel: kT e = 15 keV. 
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Table A.l. Continued. 


Star 

Instrument 

date 

TJT 

Exp. time (sec) 

S/N 

0 1 Ori C 

NARVAL 

Mar 10 th 2007 

08:28:19 

4x400 

500-1400 


NARVAL 

Mar lO* 2007 

08:59:01 

4x400 

450-1250 

.fOri A 

NARVAL 

Oct 18“ 2007 

00:52:56 

48 x 4 x 20 

780 - 990 


NARVAL 

Oct 19 th 2007 

04:35:04 

8 x 4 x 40 

1010-1080 


NARVAL 

Oct 20 lh 2007 

00:59:39 

44 x 4 x 40 

1220-1470 


NARVAL 

Oct 21 s * 2007 

23:06:15 

48 x 4 x 40 

810-1460 


NARVAL 

Oct 22”* 2007 

23:45:52 

48 x 4 x 40 

1090-1480 


NARVAL 

Oct 23 rd 2007 

23:13:45 

48 x 4 x 40 

1030-1480 


NARVAL 

Oct 24 th 2007 

23:58:34 

48 x 4 x 40 

1200-1470 

HD 57682 

ESPADONS 

May 4 th 2009 

06:05:55 

4x600 

500-1130 


ESPADONS 

May 5 th 2009 

06:22:17 

4x540 

450-1000 


ESPADONS 

May 7 th 2009 

06:46:41 

4x540 

350 700 


ESPADONS 

May 8 th 2009 

06:21:09 

4x540 

350 800 


ESPADONS 

May 9 th 2009 

06:20:36 

4x540 

300-600 

T SCO 

ESPADONS 

May 23 r0 2005 

09:14:25 

4x60 

700-1550 


ESPADONS 

Sep 20* 2005 

05:01:56 

4x30 

600 1200 


ESPADONS 

Aug 5 th 2006 

05:25:22 

4x30 

500 1150 


ESPADONS 

Mar 4* 2007 

16:23:15 

4x30 

700-1350 


ESPADONS 

Jun 29 tb 2008 

08:03:07 

4x45 

400-1050 



